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Abstract 

We establish the optimal input signaling and the capacity of MIMO channels under per-antenna 
power constraint. While admitting a linear eigenbeam structure, the optimal input is no longer diago- 
nalizable by the channel right singular vectors as with sum power constraint. We formulate the capacity 
optimization as an SDP problem and solve in closed-form the optimal input covariance as a function 
of the dual variable. We then design an efficient algorithm to find this optimal input signaling for all 
channel sizes. The proposed algorithm allows for straightforward implementation in practical systems 
in real time. Simulation results show that with equal constraint per antenna, capacity with per-antenna 
power can be close to capacity with sum power, but as the constraint becomes more skew, the two 
capacities diverge. Forcing input eigenbeams to match the channel right singular vectors achieves no 
improvement over independent signaling and can even be detrimental to capacity. 

I. Introduction 

The capacity of a MIMO wireless channel depends on transmit power constraints and the 
availability of channel state information (CSI) at the transmitter and receiver. With sum power 
constraint across all transmit antennas, the capacity and optimal signaling are well established. 
For channel state known at both the transmitter and receiver, the capacity can be obtained by 
performing singular value decomposition of the channel matrix and water-filling power allocation 
on the channel eigenvalues [1]. For Ray lei gh fading channels with state known only at the 
receiver, the ergodic capacity is obtained by sending independent signals with equal power from 
all transmit antennas [2]. 

Under the per-antenna power constraint, the MIMO capacity is less well understood. However, 
this per-antenna power constraint is more realistic in practice than sum power because each 
antenna is usually connected to a separate power amplifier on the individual RF chain. Each 
power amplifier has its own dynamic range, hence the transmitter may not be able to allocate 
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power arbitrarily among its antennas. Another appealing scenario for per-antenna constraint is 
in a distributed MIMO system, which has transmit antennas located at different physical nodes 
that cannot share power with each other. Thus understanding the capacity and optimal signaling 
schemes under per-antenna power is useful. 

The per-antenna power constraint has been investigated under different settings. In [3], the 
problem of a multi-user downlink channel is considered with per-antenna power constraint. For 
downlink broadcast, the capacity optimization problem is non-convex. It was argued that linear 
processing at both the transmitter (by multi-mode beamforming) and the receiver (by MMSE 
receive beamforming with successive interference cancellation) can achieve the capacity region. 
Using uplink-downlink duality, the boundary points of the capacity region for the downlink 
channel with per-antenna constraint can be found by solving a dual uplink problem, which 
maximizes a weighted sum rate for the uplink channel with sum power constraint across the 
users and an uncertain noise. Even though the dual uplink problem is convex, no efficient 
algorithms yet exist for solving it. 

Similar uplink-downlink duality holds with sum power constraint, in which case, the uplink 
capacity can be solved using the efficient and distributed iterative water- filling algorithm [4]. 
This algorithm is based on the closed-form solution for capacity of a single-user MIMO channel 
with sum power constraint. Yet again for per-antenna power constraint, no capacity solutions so 
far exist even for single-user systems. The lack of such a single-user solution implies full and 
centralized computation required to solve multi-user problems, including the convex problem of 
uplink capacity. 

In a previous paper [5], we have established the capacity of a single-user MISO channel (with 
single receive antenna) under the per-antenna constraint. The MISO channel admits a closed-form 
solution for the optimal input signaling as single-mode beamforming with only the beam phases 
being matched to the channel while the amplitudes are determined by the power constraint. 

In this paper, we solve the single-user MIMO capacity with per-antenna. With perfect CSI 
at both transmitter and receiver, it can be shown that channel eigen-beamforming is no longer 
optimal, as was the case with sum power constraint. We formulate the capacity optimization 
problem in the SDP framework and analyze the optimality conditions. As the optimal input is 
Gaussian, we establish in closed-form its covariance matrix as a function of the dual variable. We 
then propose a simple, iterative algorithm to find this optimal input covariance and the capacity. 
Without CSI at the transmitter (CSIT) in Rayleigh fading channels, a simple analysis shows that 
the optimal signaling scheme is to send independent Gaussian signals from all transmit antennas. 

This paper consists of 9 sections, with problem setups and preliminaries in Sections II and 
III. Sections IV- VI discuss the case of perfect CSIT with detailed solutions and algorithms. 
Section VII discusses the case of Rayleigh fading channel without CSIT. Numerical examples 
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and analysis are provided in Section VIII before conclusion in Section IX. 

For notation, we use bold face lower-case letters for vectors, bold face capital letters for 
matrices, for transpose, (.)* for conjugate, for conjugate transpose, and ^ for matrix 
inequalities (positive semi-definite relation), tr(.) for trace, diag{.} or diag(.) for forming a 
diagonal matrix with the specified elements or from the diagonal values of the specified matrix. 

II. Channel Model and Power Constraints 

A. Channel model 

Consider a frequency-flat multiple-input multiple-output (MIMO) channel with n transmit and 
m receive antennas. The channel between each transmit-receive pair is a complex, multiplicative 
factor hij. Denote the channel coefficient matrix as H of size m x n, and the transmit signal 
vector as x = [x 1 . . . x n ] T . Then the received signal vector of length m can be written as 

y = Hx + z (1) 

where z ~ CJ\f(0, 1) is a vector of additive white circularly complex Gaussian noise. Here we 
have normalized the noise power at all receivers, which can be done by absorbing the actual 
noise power into the transmit power constraint. 

Assume perfect CSI at the receiver, we will consider two cases of CSI at the transmitter: 
perfect CSIT, and no CSIT. Perfect CSIT with H known also to the transmitter can be applied 
to slow fading channels, in which channel tracking (by any mean of reciprocity or feedback) 
is possible, or to non-fading channels such as those in digital subscriber lines. No CSIT is 
applicable in fast-fading wireless channels, in which the transmitter only knows the channel 
distribution. For ergodic capacity, we model H as a random, circular complex Gaussian matrix 
with zero mean and covariance I mxn - 

The capacity of MIMO channel (1) depends on the power constraint on input signal vector x. 
In all cases, because of the Gaussian noise and known channel state at the receiver, the optimal 
input signal is Gaussian with zero mean [2]. Let Q = E[xx^] be the covariance of this Gaussian 
input, then the achievable transmission rate for a specific channel is 

r = logdet (l m + HQH f ) . (2) 

The remaining question is to establish the optimal Q that maximizes this rate or its expected 
value, hence achieves the capacity, according to the CSIT condition and a given power constraint. 

B. Power constraints 

The MIMO capacity is often studied with sum power constraint across all antennas. In this 
paper, we consider a more realistic per-antenna power constraint. For comparison, we also include 
independent multiple-access power constraint. We elaborate on each power constraint below. 
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1 ) Sum power constraint: With sum power constraint, the total transmit power from all n 
antennas is P, but this power can be shared or allocated arbitrarily among the transmit antennas. 
This constraint translates to a condition on the input covariance as 

tr(Q) < P. (3) 

This constraint allows complete cooperation among the transmit antennas. 

2) Independent multiple-access power constraint: In this case, each transmit antenna has its 
own power budget and acts independently. This constraint can model the case of distributed 
transmit antennas, such as on different wireless nodes scattered in a field, without explicit 
cooperation among them. Let Pi be the power constraint on antenna i, then the multiple- 
access constraint is equivalent to having a diagonal input covariance Q = diag{Pj}. Denote 
P = diag{Pj}, where tr(P) = P in relation to (3), then the multiple-access power constraint 
can also be expressed as 

Q ^ P. (4) 

Writing this constraint in the above semi-definite form is convenient for analyzing the capacity 
optimization problem later, but it does not alter the solution since to achieve the capacity, the 
power constraint must be met with equality and hence Q = P. 

3) Per-antenna power constraint: Here each antenna also has a separate transmit power budget 
of Pi (i — 1, . . . , n) but can cooperate with each other in terms of signaling. Such a channel 
can model a physically centralized MIMO system, in which the per-antenna constraint comes 
from the separate RF chain at each antenna. This channel can also model a distributed (but 
cooperative) MIMO system, in which each transmit antenna belongs to a sensor or ad hoc node 
distributed in a network. The distributed nodes have no ability to share or allocate power among 
themselves. The per-antenna constraint is equivalent to having the input covariance matrix Q 
with fixed diagonal values Qu = Pi. But different from the multiple-access constraint, the off- 
diagonal values of Q can be non-zero. Denote = [0 . . . 1 . . . 0] T as a vector with the i th element 
equal to 1 and the rest are 0. Then per-antenna constraint can also be written as 

ejQei<Pi, i = l...n. (5) 

This is a set of linear constraints on Q. Again expressing the constraint as an inequality instead 
of equality does not alter the optimal solution. It should be stressed that constraints on the 
diagonal values of Q are not the same as constraints on the eigenvalues of Q. 

III. ElGENBEAM STRUCTURE 

The eigen-decomposition of Q is often associated with multimode beamforming. The eigen- 
vectors of Q provide the beam directions, whereas the eigenvalues provide the power allocation 
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on these beams. Next, we will review known result for capacity with sum power, which admits 
nicely separate solutions for the beam directions and power allocation, by performing channel 
singular value decomposition. We will then discuss the beam directions for per-antenna power 
and show that channel singular vectors are no longer applicable. 

A. Review of capacities with sum and multiple access power 

With sum power constraint (3), the optimal solution is found by the well-known water-filling 
algorithm [1]. Perform the singular value decomposition of the channel as H = UhA h vJj, 
which provides well-defined singular vectors for all the non-zero singular values. Then the 
capacity-optimal covariance matrix Q* has the eigenvalue decomposition as Q* = V^AqVh, 
where the eigenvalues Aq^ are obtained through water-filling as 

A Q> j = (fi- ^H,i) + » i = l...n (6) 

for the Ah,i ^ (if Ahj = 0, the corresponding Aqj = 0). Here fx is the water level chosen 
such that tr(Q*) = P. 

Thus with sum power, the optimal solution for Q is diagonalizable by the channel right 
singular vectors. The optimal signaling is multi-mode beamforming with the beam directions 
specified by the channel right singular vectors (associated with the non-zero singular values) 
and the beam power allocation obtained by water- filling. 

With multiple-access constraint (4), the obvious solution is Q = P. In this case the optimal 
signaling is sending independent signals from different antennas, each with the constrained power. 

B. Forced beam directions with per-antenna power 

Because the channel right singular vectors provide a simple yet optimal result for the beam 
directions with sum power, we may be tempting to use the same beam directions for per-antenna 
power. But a simple analysis can show that these beam directions may not be feasible, let alone 
being optimal. That is, we may not always be able to find a Q with eigenvectors given by V H 
that satisfies the per-antenna power constraint. 

Indeed, let V H be the eigenvectors of Q for per-antenna power, then we can write Q = 
Y^i=i ^Q,i v i v l-> where are columns of V H , and Aq^ are the eigenvalues of Q. We now only 
need to find to satisfy the per-antenna power Qjj = Pj. This condition translates to 

n 

i=i 

where Vji denotes the entry of V H . Now form a new n x n matrix W with elements 
Wji = \Vji\ 2 and express the eigenvalues of Q in a vector form as Aq, then we obtain 

WA Q = p, (7) 
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where p = [P x . . . P n ] T . (Here p is a vector containing the diagonal values of P.) 

If the power constraint for each antenna is the same, that is p = ^1, the above equation has 
a unique solution of Aq = ^. This implies Q = ^I n , which is the same as the solution with 
multiple access constraint. Thus forcing the beam directions to be V H in this case is the same 
as sending independent signals from different antennas, i.e., no input optimization. 

For any other p, if W is full-rank then the eigenvalues can be found as Aq = W x p. 
The problem, however, is that the obtained Aq may be negative, thus Q may be non-positive 
semidefinite. In other words, a solution may not exist. Illustration can be found in Figure 3 of 
the numerical section for a 2 x 2 channel, in which for the infeasible cases (equation (7) does 
not admit non-negative solution), the obtained transmission rate is zero. 

IV. Capacity optimization with perfect CSIT 

In this section, we analyze the optimization problem of finding MIMO capacity with chan- 
nel known at both the transmitter and receiver. For all stated power constraints, this capacity 
optimization can be cast as follows. 

max logdet (l m + HQH f ) (8) 
s.t. <?(Q,P)<0 
Q >0, 

where g(Q, P) < refers to a power constraint as in (3), (4) or (5), and Q is Hermitian. 

Since all considered power constraints are linear in Q, the above optimization is convex with 
any power constraint. By Slater's condition [6], because of the strictly feasible value Q = P y 
which readily satisfies all power constraints, the optimal solution always exists. Thus for each 
power constraint, the problem admits a unique optimal solution for Q. For the per- antenna 
constraint (5), this optimal Q is not yet known. 

While the convex structure of this problem allows the optimal solution to be found numerically 
using available convex optimization software, such a numerical solution may be too complex for 
real-time system implementation and offers little analytical insights. Next, we will analyze this 
problem and the optimality conditions to find the optimal solution analytically. In these analyses, 
we assume that the channel H is full-rank. 

A. Capacity optimization using SDP framework 

Problem (8) can be analyzed using the SDP framework. Let D = diagjrfj} )p be a diagonal 
matrix consisting of Lagrangian multipliers d, L for the per-antenna power constraints in (5), and 
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M )p be the Lagrangian multiplier for the positive semi-definite constraint. Both D and M 
have size nxn. Then the Lagrangian for problem (8) with per-antenna power can be formed as 

£(Q, D, M) = log det (l m + HQH f ) - tr[D(Q - P)] + tr(MQ). (9) 

Note that we can form similar Lagrangian for the problem with sum power or multiple access 
constraint by replacing D with v\ n (scaled identity) or B (full matrix), respectively. These facts 
will be useful later when we compare solutions of different power constraints in Section VI-D. 
Now we will focus on per-antenna power. 

Taking the first order derivative of C in (9) with respect to Q (see [7] Appendix A.7 for 
derivatives of a function with respect to a matrix) and equating to zero, we obtain 

H f (I m + HQH f ) 1 H - D + M = 0. 

Based on the KKT conditions, we then obtain a set of optimality conditions as follows. 

H f (I m + HQH f ) 1 H = D - M (10) 

MQ = 

diagonal D >- 

Hermitian M, Q ^ 0. 

Since the problem is convex, the optimal Q is the solution to the above set of equations. Next, 
we will analyze this set of equations to first deduce a condition on the rank of optimal Q, then 
provide an equation for solving for Q. 

B. Rank of the optimal input covariance 

Following arguments similar to [5] (Appendix B), multiplying both sides of the first equation 
in (10) on the right with Q and applying the complementary slackness condition MQ = 0, we 
obtain 

DQ = H f (I m + HQH f ) 1 HQ. (11) 

At optimum, we must have D y 0. This is because each constraints in (5) must be met with 
equality, for otherwise we can always increase a diagonal value of Q and get a higher rate; hence 
the associated dual variables are strictly positive. Thus at optimum, D is full-rank, subsequently 
(11) implies that 

rank(Q) < rank(H). 

Therefore, the rank of the optimal input covariance is no more than the channel rank. In other 
words, the number of independent signal streams (or modes) should be no more than the rank 
of the channel. This result is similar to that with sum power constraint. 
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Since channel H can support at most r = min{m, n} independent modes (independent signal 
streams), the above condition implies that the rank of Q is at most r. When the rank of Q is less 
than r, that implies mode-dropping (similar to the same concept with sum power constraint). 
Since QM = 0, M is a positive semidefinite matrix in the null space of Q. The rank of M 
corresponds to the number of modes that has to be dropped for Q to be positive semidefinite. 
Suppose that the optimal solution has k modes dropped (0 < k < min{m, n}), then 

rank(M) = k, rank(Q) = min{m, n} — k. (12) 

The difference between the rank of Q and the size of Q should be stressed here. The size of 
Q is n x n. If n > m (more transmit than receive antennas), the optimal Q is inherently rank- 
deficient. In this case, even without any mode-dropping, the maximum rank of Q is m < n. 
Thus no mode-dropping does not always imply full-rank Q. Only if n < m then Q can be 
full-rank without mode-dropping. 

C. Optimality conditions with per-antenna power 

From the set of optimality conditions (10), we can obtain the following lemma. 

denote D = D _1 and define 
= HQH f 

= HDHt, (13) 
then the optimality conditions (10) imply 

(R m - F m + I m )R m = 0. (14) 

The proof is given in Appendix A. Note that both R m and F m are m x m Hermitian 
matrices and the achievable rate for each channel state in (2) can now be expressed as r = 
logdet (I m + R m ), which is a sole function of R m . 

Condition (14) provides the equation for solving for Q. To understand the meaning of this 
equation better, lets denote 

Sm Rm F m -|- I m . (15) 

Then we can show that (see Appendix A) 

HDM = S m H. (16) 

Since D is square and full rank, (16) implies rank(S m ) = rank(M) = k. From (12) and (13), 
we have rank(R m ) = rank(Q) = min{m,n} — k. From (14), we have S m R m = 0. Thus S m is 



Lemma 1. As D is full-rank and invertible, 

Rm 
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a matrix in the null space of R m , in the same way that M is in the null space of Q (but not 
necessarily spanning all the null space). In other words, R m contains the active transmission 
modes, and S m contains the modes that are dropped. Equation (14) essentially transforms the 
slackness condition from (M, Q) space to (S m , R m ) space. 
Next, we will use (14) to solve for the optimal Q. 

V. Optimal input covariance with per- antenna power 

In this section, we establish the optimal value of Q as an explicit function of the dual variable 
D, using Lemma 1. For the preliminary, let the singular value decomposition of the channel 
matrix be 



where 

• Uh is a m x m unitary matrix containing the left singular vectors, 

• Vh is a n x n unitary matrix containing the right singular vectors and 

• Eh is a m x n diagonal matrix containing the (real) singular values in decreasing order. 
It is now necessary to distinguish two cases of channel sizes: n < m and n > m. 

A. Case n < m (fewer transmit than receive antennas) 

For n < m, the optimal Q solution can have up to all n modes. In Lemma 1, however, F m 
is not full-rank. To proceed, we need to convert this matrix to full-rank as follows. 

In (17), for n < m, we can write S H = [S„ 0„ i?n „„] T , where S n is a n x n diagonal matrix 
containing the real (non-zero) singular values of H. Now define 



H = U H S H V 



H 



(17) 



K = V H S„V : 



(18) 



then K is square and full-rank. From Lemma 1, we can derive the following result. 



Lemma 2. Define two n x n matrices as 



R„ = KQK f 
F n = KDK f . 



(19) 



then for n < m, the optimality conditions (10) imply 



(Rn — F n + I„)R n — 0. 



(20) 



We get an equation similar to (14), but here for n < m, F n is full-rank. 
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Proof: Multiplying (14) on the left with W and on the right with H, and noting that 
HtH = KKt, we obtain 

K[(R n -F n + I n )R„] Kt = 0. 

Since for n < m, K is square, full-rank and hence is invertible, the above equation is equivalent 
to (20). ■ 
We now analyze equation (20). This equation can be written as R^+R„ = F n R n . This equality 
implies that F n R„ is Hermitian and has the same eigenvalue decomposition as R£ + R ra , which 
has the same eigenvectors as those of R n . This is possible only if R n and F n share the same 
eigenvectors for the non-zero eigenvalues. Now since for n < m, F n is full rank and Hermitian, 
it has a unique eigenvalue decomposition (up to any multiplicity of eigenvalues). From (19), 
rank(R„) = rank(Q) = n — k. Equation (20) then implies that R„ must span n — k eigenspaces 
of F n : specifically, the n — k eigenvectors corresponding to the non-zero eigenvalues of R„ are 
the same as n — k eigenvectors of F n . Define 

S n = R„ — F n + I n , (21) 

then equivalently, (20) implies that S n spans the other k eigenspaces of F n : the k eigenvectors 
with non-zero eigenvalues of S n are the same as the other k eigenvectors of F n . 

Intuitively, this result can be interpreted as follows. Note that we can write F n — I n = R n — S n . 
The matrix F n — I n may contain some positive and some non-positive eigenvalues. Then R n is 
the portion that contains only the positive eigenmodes, and (— S„) is the portion that contains 
only the non-positive eigenmodes. As such, both R n and S„ are positive semidefinite matrices 
and are orthogonal to each other. Here R„ contains the (positive) n — k transmission modes, 
while S n contains the k modes that are dropped. 

Based on this analysis, we can obtain the optimal value for Q as a function of D as follows. 

Theorem 1. For n < m, K as defined in (18) is full-rank and invertible. Denote K = K _1 , 
then for a given D >- 0, the optimal Q satisfying the optimality conditions (10) is given by 

Q* = D - KK f + Z, (22) 

where Z = KS n K^, and (— S„) is obtained as the non-positive eigenmodes o/KDK* — l n . 

Proof: From (19), we have Q = KR n Rt. Result (22) then follows directly from (20), (21) 
and the associated analysis. ■ 
Theorem 1 gives the solution for Q in terms of the dual variable D. Here K is a function 
of the channel as defined in (18), while S„ is determined from F n which is a function of D as 
in (19). Note that since R„ contains only the positive eigenmodes of F n — I n , the optimal Q* 
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as formed in (22) is always positive semidefinite. Thus the only step left is to find the optimal 
dual variable D such that Q* satisfies the power constraint of diag(Q*) = P. 

To find the optimal D, at this point, we need to use an iterative algorithm which we will 
discuss in Section VI. 

B. Case n > m (more transmit than receive antennas) 

For n > m, the optimal Q is inherently rank-deficient since the channel can support at most 
m modes, which is smaller than the number of transmit antennas. For this case, we need to 
further decompose channel H as follows. 

In (17), for n > m, we can write S H = [S m O m ,n-m], where S m isamxm diagonal matrix 
containing the singular values of H. Now separate the right singular vectors of H as 

V H = [Vi V 2 ], Vi = first m columns ^ 
V 2 = last n — m columns. 

Here Vi contains the basis for the row space of H, while V 2 contains the basis for the null 
space of H. Only Vi is unique, but V 2 can be any basis matrix spanning the null space of H 
(i.e. of Vi). Note that we can also write H = UHS m V{. Now lets denote the channel "inverse" 
H as 

H = y 1 H m - 1 Ul i , (24) 

then HH = I m . The use of H will become apparent later. 

Next, from (10), multiplying both sides of the first equation on the left with ~V\ and on the 
right with Q, then applying the second equation, we get 

V^DQ = 0. (25) 

Equation (25) places a constraint on the rank of Q as a direct consequence of m < n. 

In this case, F m as defined in (13) is full rank and Hermitian, hence it has unique eigenvalue 
decomposition (up to any multiplicity of eigenvalues). Equation (14) of Lemma 1 then implies 
R m and S m share eigenvectors with F m , where R m spans m—k eigenspaces of F m , and S m spans 
the other k eigenspaces of F m . Specifically, R„ contains the positive eigenmodes of F m — I m 
while (— S m ) contains the non-positive eigenmodes. (The negation is just for convenience so 
that S m is positive semidefinite.) 

Now from (13) and (15), we obtain 

H(Q - D)Ht = S m - I m . (26) 

Different from (20) of Lemma 2, the above equation is under-determined for Q. In order to find 
the unique optimal Q*, we need to combine (26) with the rank condition (25). Based on these 
two equations, we can obtain the optimal Q* as follows. 
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Theorem 2. For n > m, establish Vi, V 2 as in (23) and H as in (24). Then for a given D >- 0, 
?/ze optimal Q satisfying the optimality conditions (10) z's gz'ven /Vy 

Q* = D - HH t + Z - X, (27) 

where Z = HS m IT, and (— S m ) zs obtained as the non-negative eigenmodes o/HDtf — I m . 
//ere X is a Hermitian matrix given as 

X = V 2 AV^ + VxB V\ + V 2 B t Vj , (28) 

where A is a (n — m) x (n — m) Hermitian matrix and B is a m x (n — m) matrix given by 

B = Vj (Z — HH^ D y 2 (V^DV^ 1 , 

A = (l n _ m - BtvlDV 2 ) (vjDV 2 ) _1 . (29) 
The proof is in Appendix B. 

Equations (27)-(29) of Theorem 2 give the solution for Q in terms of the dual variable D in 
a form similar to (22) of Theorem 1. However, there are extra terms here involving V 2 which 
spans the null space of the channel. In (27) and (29), Z is a function of S m which is determined 
from F m , which in turn is a function of D. Thus again the only remaining step is to find a 
diagonal D >- such that diag(Q*) = P. 

C. The duality gap 

To better understand the results of Theorems 1 and 2, we now analyze the duality gap. The 
solutions of Q* in (22) and (27) always satisfy MQ = and Q )p as depicted in the set 
of optimality conditions (10). For a given dual variable D, such a Q* is precisely the optimal 
primal variable to establish the Lagrange dual function from the Lagrangian (9) as 

£*(D) = max£(Q,D,M) 

V ' Q,M V ' 

= logdet (I n + HQ*H f ) - tr[D(Q* - P)]. 

The optimal dual variable D* is then the value that minimizes the above dual function. Solving 
for D* by directly minimizing £*(D), however, is not simple because the derivative of Q* with 
respect to D is complicated and may not even be derivable in closed-form. However, since 
problem (8) is convex and satisfies Slater's condition, it has zero duality gap at optimum. For a 
given D, the duality gap is 

G(D) = -tr[D(Q* - P)]. (30) 
Any algorithm that has this duality gap approach zero will converge to the optimal value. 
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VI. Algorithm for finding the optimal Q 

In this section, we investigate the remaining question of finding a diagonal matrix D y such 
that the solutions in Theorems 1 and 2 satisfy diag(Q*) = P. There appears to be no closed-form 
analytical solution for such a D. Fortunately, equations (22) and (27) in these theorems suggest 
a way to compute D iteratively. In the follows, we design an iterative algorithm for finding the 
optimal D, and hence optimal Q*, for each case of n < m and n > m. Then we integrate both 
cases in a main program for any channel size. 

A. Iterative algorithm for finding Q* when n < m 

Equation (22) in Theorem 1 suggests a simple iterative algorithm for finding the optimal D. 

First we need to choose an initial point D . This point can be chosen arbitrarily as long as 
it satisfies D >- 0. For potentially faster convergence, we follow the mode-dropping approach 
similar to water filling. For the initial point, we assume that there is no mode-dropping. Denote 
G = KK+, then based on (22), we can just simply choose diagonal matrix D as 

D = P + diag(G). (31) 

This solution always satisfies D y since G as a positive semidefinite matrix has non-negative 
diagonal values. At this step, we can perform a quick check to see if D — G is positive 
semidefinite. If it is, then this value is the optimal input covariance, i.e. Q* = D — G, and no 
iteration is needed. 

If (D — G) is non-positive semidefinite, then we adjust D using an iterative procedure as 
follows. Say at iteration i (i > 0), we have obtained Dj. Then we can form F Hii , S Hii and Qj as 

F n>i = KD.Kt 
— S n: i = non-positive eigenmodes of (F n i — I n ) 

Qi = Dj - G + Zi. (32) 

The Qi as computed in (32) is always positive semidefinite (as a consequence of Theorem 1). 
The term Z; > can be thought of as the adjustment at step i to make Qj ^ 0. But it does not 
guarantee that the diagonal of Qj will be P. From (32), noting that Dj is diagonal, we update 
Dj + i by the difference between the diagonal of Qi and P as 

D m = D i + P-diag(Q,). (33) 

Iteration stops when the diagonal values of Qj is close to P within an acceptable tolerance. In 
implementation, we choose to stop when the duality gap (30) satisfies |tr[D(Q — P)]| < e. Since 
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Algorithm 1 drop-rank-n(n, D , K, K, G, P, e): Iterative search for Q* when n < m. 
Require: D diagonal y 0, P diagonal y 0, e > 0. Also K = K 1 and G = KK 1 ". 

l: Initialize % — 0. (iteration count) 

2: Initialize A = 1 + e. (loop terminating variable) 

3: while (A > e) do 

4: Form Fj = KDjK^ — I n . (note that — I n is included here) 

5: Perform the eigenvalue decomposition F, = UfA f U f . Let k be the number of non- 
positive eigenvalues of Fj. 
6: Form Sj = -UpApUp* where 

Ap is the k x k diagonal matrix of all k non-positive eigenvalues of Fj, 
U F consists of the corresponding k eigenvectors. 
7: Form Zi = KS^K^ 
8: Form Qj = Dj - G + 
9: Form D, +1 = + P - diag(Q,). 
10: Compute A = |tr[Dj(Qj - P)]|. 
11: + 

12: end while 

13: return Dj and Qj. 



problem (8) is convex and satisfies Slater's condition, this stopping criterion always guarantees 
the optimal solution. The iterative procedure for finding Q* when n < m is summarized in 
Algorithm 1, drop-rank-n(-). 

B. Iterative algorithm for finding Q* when n > m 

For n > m, we rely on equation (27) in Theorem 2 to design a similar algorithm. Let 
G = HH^, (note that we reuse the symbol G here but there should no confusion based on the 
channel size). Again we can start with the initial value D in (31) or with any arbitrary diagonal 
Do y 0. In this case, however, iteration is always necessary (except in the unlikely event that 
the algorithm starts with the optimal D). 
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Algorithm 2 dwp-rank-m(m, D , H, H, G, Vi, V 2 , P, e): Iterative search for Q* when n > m. 
Require: D diagonal >- 0, P diagonal y 0, e > 0. Also Vi and V 2 are related to H as in 
(23), H is defined in (24) and G = HHt 
Initialize % — 0. (iteration count) 
Initialize A = 1 + e. (loop terminating variable) 
while (A > e) do 

Form Fj = HDjtP — I m . (note that — I m is included here) 
Perform the eigenvalue decomposition Fj = UfA f U f . 
Form Sj = — UpApUp^ that contains all k non-positive eigenmodes of Fj. 
Form Zj = HSjHt. 

Form Dj = diag{(Ajj)~ 1 }, j = l...n. 
Form B, = Vj (Z, - G) D,V 2 (v^DjV 2 ) _1 . 

Form A, = (l n _ m - BjvlDjV 2 ) (v^D,V 2 ) _1 . 
Form X, = V 2 A l Vt + V 1 B l \t + V 2 BtvJ. 
Form Qj = Dj - G + Zj - Xj. 
Form D J+1 = Dj + P - diag(Q,). 
Compute A = |tr[Dj(Qj - P)]|. 

end while 
return Dj and Qj. 



At step i (i > 0), having obtained Dj, we compute 

F m ,j = HDjH t 
— S m>i = non-positive eigenmodes of (F m>i — I m ) 
Zj = HS mj jH^ 

X, = V 2 A i V^ + V^Vt + V 2 BjvI 

Qi = Dj - G + Zj - Xj, (34) 

where Aj and Bj are computed from Dj and Zj as in (29). The Qj as computed in (34) is again 
always positive semidefinite and has rank at most m. Then we update D i+1 in the same way as 
in (33), and stop when the duality gap |tr[D(Q — P)]| is sufficiently small. This procedure is 
summarized in Algorithm 2, drop-rank-m(-) . 
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Algorithm 3 opt-cov(H, {Pi}, e): The main program to find Q* for a given channel H. 
Require: Channel H full-rank, per-antenna power P, L real and Pj > 0, % — 1 . . . n, precision 

e > 0. 

l: Let (to, n) = size(H), 

to = number of rows or receive antennas, 

n = number of columns or transmit antennas. 
2: Perform the singular value decomposition H = UhS h V^, where 

U H : m x to unitary matrix, 

V H : n x n unitary matrix, 

S H : to x n diagonal matrix with diagonal entries as real, non-zero singular values cr H ,i 
in decreasing order. 

3: if n > to then 

4: Form V 2 from the last n — m columns of V H . 

5: Form Vi from the first to columns of Vh- 

6: Form S m = diagjo-j^}, to x to diagonal matrix. 

7: Form H = ViS m U H . 

8: Form G = ViE^Vj. 

9: else 

10: Form S n = diag{o- H ,i} 5 ^ x n diagonal matrix. 
11: Form K = V H S n V^. 

12: Form S n = diagjoj^}, n x n diagonal matrix. 
13: Form K = V H S n V^. 
14: Form G = V H S^V^. 
15: end if 

16: Form P = diag{P}, n x n diagonal matrix. 
17: Form D = P + diag(G). 
18: if n < to then 

19: (D, Q) = drop-rank-n(n, D , K, K, G, P, e) 
20: else 

21: (D, Q) = drop-rank-m(m, D , H, H, G, Vi, V 2 , P, e) 
22: end if 
23: return Q. 



The integrated main program for finding the optimal Q* for all channel sizes is summarized 
in Algorithm 3, opt-cov(-). 
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Fig. 1. Typical convergence for a 3 x 3 channel at SNR = -5dB and e = 10 



C. Convergence analysis 

1) For n < m: The following Corollary shows convergence of Algorithm 1. 

Corollary 1. In Algorithm 1, Dj is decreasing in i. This algorithm always converges to the 
optimal point. 

Proof: Consider the iterative steps (32) and (33) with initial condition (31). We will show 
that Dj is decreasing in i by induction as follows. 

First, with i = 0, from (31) and (32) we have Q = P + Z + diag(G) - G. Thus from (33), 

Di = D — diag(Z ). 

Since Z ^ 0, then diag(Z ) ^ 0, and thus Di ^ D . 

Now assume that t)j+\ =4 D? for some j > 0. From (32), we then have F nj - + i ^ F n j. Since 
F n — I n = R n — S„, this implies S nj — S n> j+i =4 R„,j — R nj+ i. Multiplying both sides on the 
left with S n j and noting that S„jR nj = 0, we obtain S n j (S nj — S nj -+i) ^ — SnjRnj+i ^ 0. 
Since S n j )p 0, this implies S n j ^ S„ Now multiplying both sides on the left with K and 
on the right with K^, we then have Zj =4 Zj+i- 

From (33) and (32), we can write 

D j+1 = P + diag(G)-Z j . 

Thus with Zj =4 Zj+i, then t)j+2 =4 This completes the induction. 

The sequence Dj is decreasing in i and is lower bounded as Dj y 0. Thus D ?; must converge. 
Based on (33), then Qj must also converge to an optimal value with diag(Q) = P. In other 
words, the algorithm always converges to the correct optimal value. ■ 
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Iteration Iteration 
Fig. 2. Typical convergence for a 2 x 4 channel at SNR = -5dB, e — 1CP 6 and random starting point Do- 

Note that convergence holds for any starting point, not just the value of Do in (31). Figure 
1 shows a typical convergence behavior of Dj and diag(Qj) as well as the objective function 
value and the duality gap for n = m = 3 (the channel is generated randomly according to the 
circularly complex Gaussian distribution). The convergence rate appears to be exponentially fast. 

2) For n > m: Now consider Algorithm 2 with iterative steps (34) and (33) and initial 
condition (31). The algorithm in this case is more complex. Here Dj is not always decreasing in 
i, but some of its diagonal elements are decreasing and others increasing. However, we observe 
that D j7r is decreasing in i for some diagonal matrix ir with only 1 and —1 on the diagonal such 
that for the initial step, diag(Q )7r )p Ptt (or after a number of iterations K, diag(Qx)vr >p Pit). 
That is Dj7r )p D i+1 7r. The detailed analysis is given in Appendix C. 

Given Dj7r is decreasing in i, then Dj converges, which implies from (33) that diag(Qj) 
converges to P. The convergence also holds for any starting point D . 

Figure 2 shows a typical convergence behavior for a random channel with n = 4 and m = 2 
at low SNR (— 5dB). Again the convergence rate appears to be exponentially fast. Numerical 
simulations show that higher SNR generally leads to faster convergence. 

D. Connections with water-filling and multiple-access 

As noted in Section IV-A, all 3 power constraints can be analyzed in the SDP framework (8). 
The only difference is in the dual variable associated with each power constraint in the Lagrangian 
(9). This dual variable for the sum, per-antenna and multiple-access constraint respectively is 
equal to v\ n (scaled identity), D (square diagonal) and B (square full matrix). This change of 
dual variable, however, makes it possible to find closed-form solutions for the sum and multiple 
access constraints. Next we will discuss these solutions and the implication on input eigenbeams. 



SUBMITTED 



19 



1 ) Optimal solutions for sum and multiple access constraints: For the sum power constraint, 
all results of Theorems 1 and 2 apply by replacing D with v\ n . This change, however, makes 
it possible to find closed-form solutions for v even with mode-dropping. In particular, with sum 
power constraint, the result of Theorem 1 becomes 



where (— S n ) now contains the non-positive eigenmodes of z/ _1 KK^ — I n . Equation (35) then 
implies that Q has the same eigenvectors as those of K, which are V H (the right singular vectors 
of H). Similar eigenvector analysis holds for Theorem 2. These eigenvectors then lead directly 
to the water-filling solution (6) for the eigenvalues of Q, where p = I /v. As such, a closed-form 
solution for the optimal dual variable v* can be obtained as 



Similarly, for the multiple access constraint, results of Theorems 1 and 2 also hold with D 
being replaced by a full nxn positive semidefinite matrix B. In this case, the optimal co variance 
is Q = P and there is no mode dropping (S„ = and S m = 0). Thus for n < m, we can easily 
identify the optimal dual variable as B* = (P + G) _1 . For n > m, closed-form expression for 
B* can also be derived from (26)-(29). 

2) Comments on input eigenbeams: In both cases of sum and multiple access constraint, 
the eigenvectors of the optimal Q are independent of the dual variable. For the per-antenna 
constraint, however, eigenvectors of Q depend on the dual D as shown in Theorems 1 and 2. As 
a consequence, the optimal beamforming directions with per-antenna constraint are not the right 
singular vectors of the channel (as discussed in Section III). For MISO channels, the rank-one 
optimal Q is shown to have its eigenvector matched to the phase but not the amplitude of the 
channel vector [5]. For MIMO channels, the relation is more complicated as the eigenvectors of 
Q depend on multiple factors: the channel H, the per-antenna constraint P and also the SNR. It 
does not appear feasible to find the optimal eigenbeams and power allocation separately. Instead, 
the proposed algorithms establish the optimal eigenbeams and power allocation together in the 
resulted optimal Q. 

VII. Fading MIMO capacities with no CSIT 

In this section, we consider the case of fading channel with no CSIT. The transmitter does 
not know the channel realization H but only knows its distribution as a circularly complex 



Q = v~ l l n - KKt + KSnKt, 



(35) 




where K < min (m, n) is the number of active modes (assuming Ah,« are in decreasing order). 
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Gaussian matrix with zero mean and covariance l nxm . In this case, we need to consider the 
ergodic capacity. For all power constraints, capacity optimization can now be cast as follows. 

max E u [logdet (l m + HQH f )] (36) 

s.t. g(Q,P)<0, Q>0, 

where g(Q, P) < refers to a specific power constraint as in (3), (4) or (5). Different from (8), 
the optimal Q is no longer a function of H but only of its distribution. 

With sum power, the capacity of MIMO fading channel was established in the seminal paper 
by Telatar [2]. The optimal covariance of the Gaussian transmit signal is Q = -^1, implying that 
each antenna sends independent signal with equal power. With multiple-access constraint, the 
transmit covariance again has to be Q = P. 

We have established the ergodic capacity with per- antenna power for the MISO channel earlier 
in [5] and showed that it is the same as capacity under multiple-access constraint with optimal 
Q = P. The same conclusion, in fact, holds for MIMO channels. For the sake of completeness, 
we provide the analysis with a shorter proof for MIMO channels here. We use an approach similar 
to [2] but only apply the following set of diagonal matrices instead of all unitary matrices. Define 
n^) as an n x n diagonal matrix with all 1 on the diagonal except a (—1) at location j. Then 
HnW has the same distribution as H and hence 

E H [logdet (1 + HQH f )] = E u logdet (l + ffll^Qn^H 1 ") . 

Here the set of 11^ matrices preserves the diagonal values of H^QIl^ T . We can then apply 
the following inequality based on the concavity of log det function: 



r = ~ Yl Eu logdet ( x + Hn^Qn^H 1 ") 

1« .g det ( 1 + H ( - II ' QH ' ' 



< E 



Ht 



= E n [log det (l + HPH f )]. 

Hence the capacity is achieved by setting Q = P. Since the set of diagonal matrices used 
here is a subset of all unitary matrices used in [2], it follows immediately that ergodic capacity 
with per-antenna power is smaller than or equal to ergodic capacity with sum power. If and only 
if the per-antenna constraint is P = — I n , then the two capacities are equal. 

VIII. Numerical Results and Analysis 

In this section, we provide numerical examples to illustrate MIMO capacity with per-antenna 
power constraint. We will only show results for perfect CSIT, in which per-antenna power has 
significant impact. 
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Fig. 3. Capacities of a 2 x 2 channel with perfect CSIT at SNR = OdB with P = diagfA, 1 - Pi}. 



A. Impact of per-antenna power for a given channel 

First, to see the impact of separate power constraint on each antenna, we vary the constraint 
P = diag{Pi, P 2 } for a 2 x 2 channel such that the total power P 1 + P 2 = 1 (the SNR = 
OdB). Figure 3 shows a typical capacity plot for a channel randomly generated according to the 
circularly complex Gaussian distribution (as in Rayleigh fading). The specific channel for this 
figure is 

0.0541 - 0.4066i -0.4339 + 0.0033i 
-1.3200- 0.1872i 0.8269 - 0.0279i 
The figure shows that input power constraint can affect the capacity significantly. For example, 
even at equal constraints Pi = P 2 = 0.5, the capacity with per antenna power is significantly 
higher than with multiple access constraint, but is still lower than with sum power. Sum power 
constraint, however, does not result in equal power at each antenna at this point. 

Also plotted is the achievable rate for per-antenna constraint using forced eigenbeams as the 
channel right singular vectors as in (7). The 3 lines for capacity with sum power, capacity with 
per-antenna power and rate achievable with forced eigenbeams meet at a single point when the 
per-antenna constraint is such that it coincides with the optimal power per antenna under sum 
constraint. Other than this point, the rate with forced eigenbeams is always smaller than capacity 
with per-antenna power and goes to zero at certain P because no solutions exist. 

B. Average impact: equal or unequal constraint? 

Next, to see the average impact of each power constraint, we plot the ergodic capacity for 
Rayleigh fading channels with perfect CSIT. Figure 4 shows the ergodic capacity of a 4 x 2 
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Fig. 4. Ergodic capacities of a 4 x 2 Rayleigh fading channel with perfect CSIT and P = ^In- Only for this equal constraint 
P that the rate with forced eigenbeams is as large as the capacity with multiple access constraint (i.e. no input optimization). 

channel, obtained by averaging over 1000 random channel realizations at each SNR. In this figure, 
the per-antenna power has equal constraint at each antenna, i.e. P = —I n . We can observe that 
the capacity with per-antenna power is then almost as high as that with sum power. However, 
the rate with forced eigenbeams is exactly the same as capacity with multiple access constraint 
(no input optimization), which agrees with the analysis in Section III. 

If we make the power constraint at each antenna unequal, the difference in capacity with 
per-antenna power and with sum power is more pronounced, as seen in Figure 5 for a 3 x 3 
channel with power at antenna k proportional to k 2 . The more skew the constraint, the more 
different the ergodic capacity becomes under per-antenna and sum power constraint. 

Figure 5 also shows the rate with forced eigenbeams at almost zero for all SNR. Even though 
this plot emphasizes the impact by skewing the per-antenna constraint significantly, simulations 
show that even if the constraint is only slightly different from — I n , rate with forced eigenbeams 
falls strictly below rate with no input optimization. This relation readily follows from (7), since 
for P 7^ — I n , the probability that equation (7) has no non-negative solution is strictly non-zero. 
As the constraint becomes more skew, the average rate with forced eigenbeams approaches zero. 

C. Impact of channel size 

Simulations also show that per-antenna power has higher impact on ergodic capacity when 
n > m. If the power at each antenna is constrained to be the same (P = then for n < m, 

the ergodic capacities under the 3 power constraints differ only slightly at low SNR and approach 
each other as the SNR increases. For n > m, however, the ergodic capacity with per-antenna 
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Fig. 5. Ergodic capacities of a 3 x 3 Rayleigh fading channel with per-antenna constraint at antenna k proportional to k 2 . 
The capacity with per-antenna power approaches capacity with multiple access at high SNR only for n < m, but remains 
significantly higher at all SNR for n > m. The rate with forced eigenbeam is almost zero for all channel sizes. 

power remains significantly higher than that with multiple access and is relatively close to the 
capacity with sum power for all SNR. 

On the other hand, if the power at each antenna is constrained differently from each other, then 
capacity with per-antenna power can be significantly different from capacity with sum power for 
all channel sizes. As the SNR increases, capacity with per-antenna power approaches capacity 
with multiple access constraint for n < m, but remains higher for n > m. 

Last, we have also verified that for the case of single receive antenna (m = 1), the proposed 
algorithm converges to the same result as the closed-form solution in [5]. 

IX. Conclusion 

We have established the MIMO capacity under per-antenna power constraint with perfect 
transmitter channel state information. The optimal input covariance matrix is no longer diag- 
onalizable by the channel right singular vectors as with sum power constraint. We solve in 
closed-form the optimal covariance matrix as a function of the dual variable. We then design an 
efficient algorithm to find this dual variable and hence the optimal input signaling. Simulation 
results show that per-antenna constraint can impact the capacity significantly. The impact is 
more pronounced in channels with more transmit than receive antennas. In all cases, to maintain 
realistic per-antenna power requirements, solutions using channel eigenbeams produce rate equal 
or worse than that using independent signaling (no input optimization). Optimal signaling can 
achieve a capacity close to the idealistic capacity with sum constraint and significantly higher 
than capacity with independent signaling. 
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Appendix 

A. Proof of Lemma 1 

In this proof, we make repeated use of the following identity expansion: 

I m = (I m + HQHt) (I m + HQHt) 1 . (37) 

Note that the order of the two factors in this expansion is interchangeable. 

To prove (14), multiplying the first equation in (10) on the right with QHt and on the left 
with HD, and noting that MQ = 0, we get 

HDH f (I m + HQH f ) 1 HQH f = HQH f . 

Now subtracting both sides by HDHt, then applying identity expansion (37), this equation 
simplifies to 

— HDH t (l m + HQH f ) 1 = HQH f - HDHt. 

Next adding both sides with I m and again using identity expansion (37), we get 

(I m + HQH f - HDH f ) (I m + HQH f ) 1 = I m + HQH f - HDHl 

Denote S m = I rn + HQHt — HDH t similar to (15), the above equation becomes 

S m (l m - (l m + HQHt)~ 1 ) =0. 

Then applying identity expansion (37) once more, we obtain 

S m (HQHt) (I m + HQHt)- 1 = 0. 

Since (l m + HQHt) 1 j s a full-rank square matrix, the above equation is equivalent to (14). 
To show (16), from the first equation in (10), subtracting HtH from both sides, we have 

Ht (I m + HQHt) 1 _ Im j H = D - M - HtH. 
Using identity expansion (37), we obtain 

H f (I m + HQH f ) 1 HQH f H = H f H - D + M. 
Now replacing a part of the left expression by the first equation in (10), we have 

(D - M)QHtH = HtH - D + M. 

But MQ = 0, hence we have 

M = DQH f H + D - H f H. 
Multiplying both sides on the left by HD, we obtain (16). 
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B. Proof of Theorem 2 

From equation (26), multiplying on the left with H and on the right with H^, we get 

VM(Q - D)V 1 Vl = HS m Ht - Ffflt 

Denote Z = HS m H^ and G = HHl Now multiply the above equation on the left with V{ and 
on the right with Vi, and noting that VjGVx = S H 2 , we get 

vI(d + z-q)v 1 = Sh 2 - 

Since V 2 Vi = 0, the above equation implies that 



D + Z - Q = [Vi V 2 ] 



"s H 2 B " 




"vi" 


Bt A 




.vj_ 



= V^M + V 2 AVl + ViBVj + V 2 B^V{ 

for some Hermitian (n — m) x (n — m) matrix A and some m x (n — m) matrix B. Thus we 
can write Q as in (27). 

The remaining question is to find A and B such that the rank condition in (25) is satisfied. 
To do this, multiplying equation (27) on the right with DV 2 and on left with either \\ or V|, 
and noting that V 2 Z = V 2 G = 0, we obtain respectively 

= I n _ m - AV|DV 2 - BtvjDV 2 , 
= VJ(Z - G)DV 2 - BV|DV 2 . 

Noting that V2DV 2 is full rank and invertible, we then obtain (29). 

C. Convergence analysis for n> m 

In this Appendix, we analyze Algorithm 2 for n > m to support the observation that Dj7r 
decreasing in i, where n is a diagonal matrix with only 1 and —1 on the diagonal such that 
diag(Qo)vr fc= Pvt. 

To show Di+i7r ^ Dj7r, we use induction. For % — 0, from (33), we have 

Eh =D + P-diag(Q ). 

Since diag(Q )7r ^ Pit, it then follows that T>i-k ^ D vr. 

Now suppose that the inequality holds for % = j, such that Dj+i7r ^ Dj7r. From (34), we 
then have 

D i ~ = (Qj ~ Qj'+i) ~ ( Z i _ Z i+i) + ( X i ~ X i+i) ■ 
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First multiplying both sides on the right with n and applying Dj + i7r Dj7r, we obtain 

(Zj - Z j+ i) 7T - (Xj - Xj + i) 71 ^ (Qj - Q j+ i) 7T. 

Now multiplying both sides on the left with V 2 V 2 Dj )p and noting that V^D^Qj = 0, we get 
V 2 V 2 D, [(Zj - Z i+1 ) vr - (X, - X i+1 ) tt] ^ -V 2 V 2 D,Q, +1 7r. (38) 
Since V 2 V 2 Dj ^ 0, this can lead to 

(Zj - Z j+1 ) tt - (X, - X j+1 ) vr ^ -Q i + 1 7r. (39) 

Note that since V 2 V 2 Dj is rank deficient, strictly (38) need not to always imply (39). However, 
we find numerically that the premise of Dj7r decreasing always holds either from the initial 
iteration or after a finite number of iterations, regardless of the initial point D . If it holds after 
a number of initial iterations, we can reset Q vr. For this reason, we conjecture that (39) always 
holds true after some iterations. 

Next multiplying both sides of (39) on the left with V 2 V|D J+1 )? and noting that V"|Dj +1 Qj +1 = 
0, we get 

V 2 V|D J+1 [(Zj - Z j+1 ) ix - (X, - X j+1 ) vr] ^ 0, 
which similarly leads to 

(Zj - Z j+1 ) tt - (Xj - X j+1 ) tt ^ 0. (40) 
Now from (33) and (34), we have 

D m = P + diag(G) - diag (Z, - X) . (41) 
Combined with (40), since diagonal of a positive semidefinite matrix is non-negative, we have 

D j+2 7r - D i+ i7r = diag (Zj - Xj) vr - diag (Z j+1 - X j+1 ) tt ^ 0, 
which completes the induction. 
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